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We introduce a semistochastic implementation of the power method to compute the dominant 
eigenvalue and eigenvector of very large matrices. The method is semistochastic in that the matrix 
multiplication is partially implemented numerically exactly and partially in the sense of expectation 
value only. Compared to a fully stochastic method, the semistochastic approach significantly reduces 
the computational time required to obtain the eigenvalue to a specified statistical uncertainty. This 
is demonstrated by the application of the semistochastic quantum Monte Carlo method to systems 
with a sign problem, the fermion Hubbard model and the carbon dimer. 

PACS numbers: 02.70.Ss, 31.10.+Z, 31.15.V-, 71.15.-m 
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Introduction. — Consider the computation of the dom- 
inant eigenvalue of a matrix of order N so large that a 
general N x N matrix cannot be stored in full. Trans- 
formation methods cannot be used in this case, but one 
can still proceed with the power method, a.k.a. the pro- 
jection method, by virtue of the fact that no more is 
required than the result of multiplication of an arbitrary 
vector by the matrix. For even larger N, when even a sin- 
gle vector can no longer be stored, Monte Carlo methods 
can be used to represent stochastically both the vector 
and multiplication by the matrix. This operation suffices 
to implement the power method to compute the domi- 
nant eigenvalue and averages involving its corresponding 
eigenvector. 

In this paper, we propose a hybrid method consist- 
ing of numerically exact representation and multiplica- 
tion in a small deterministic subspace, complemented by 
stochastic treatment of the rest of the space. This semis- 
tochastic projection method combines the advantages of 
the deterministic and the stochastic projection methods 
- it greatly reduces the statistical uncertainty of averages 
relative to stochastic projection while allowing N to be 
large. These advantages are realized if one succeeds in 
choosing a deterministic subspace that carries a substan- 
tial fraction of the total spectral weight of the dominant 
cigenstate. 

Semistochastic projection has numerous potential ap- 
plications: transfer matrix Monte Carlo [lj for classi- 
cal statistical mechanical systems, quantum Monte Carlo 
(QMC) 2|-|j], and the calculation of subdominant eigen- 
values [5|. 

In this work we apply the semistochastic method to 
compute the ground state energy of quantum mechanical 
Hamiltonians represented in a discrete basis. In this con- 
text, deterministic projection is known as Full Configura- 



tion Interaction to chemists and as exact diagonalization 
to physicists, whereas stochastic projection is the essence 
of various projector QMC methods 0, Q ■ Hence, semis- 
tochastic projection shall be referred to as SQMC. The 
benefit of SQMC over the corresponding QMC method is 
large in many systems of interest since the Hartrce-Fock 
determinant, augmented by a small set of additional de- 
terminants, represents a significant fraction of the total 
spectral weight of the ground state wavefunction. 

The Hamiltonians for the systems considered here suf- 
fer from a sign problem, i.e., no sign changes of basis 
states can be found that render all elements of the matrix 
and the desired eigenstate non-negative. Until recently, 
projector QMC had been used used only for systems that 
do not have a sign problem 0, 0], or with an uncon- 
trolled, variational fixed- node approximation [(ij]. The 
recent breakthroughs of Alavi and coworkers with their 
FCIQMC method 0] and its initiator extension || , have 
enabled the treatment with a controllable bias of matri- 
ces with a sign problem. Consequently, the stochastic 
method to which we compare SQMC is essentially the 
same as the initiator FCIQMC (i-FCIQMC) method of 
Alavi, but with some minor differences explained in the 
text. 

Theory. — We start from an N x N Hcrmitian matrix 
H, with eigenvalues Eq < E\ < ■ ■ ■ < E^-i- In our case, 
H is the Hamiltonian represented in an orthonormal basis 
{|0i },..., \4>n)}- To obtain the lowest eigenvalue E , and 
its eigenvector ipW with components vbf^ = (<^,-|^>( )), we 
first invert, shift and scale the Hamiltonian matrix: 

P=1 + t(E t 1-H), (1) 

where Et is a running estimate of Eq. 

For the case Et = Eq, P has unit eigenvalue. If 
r < 2/(-E;v— l — Eq), then the unit eigenvalue is the 
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dominant one, and the E+-depcndence in P ensures 
that the power method iterates are bounded and non- 
vanishing. When multiplication by P is performed deter- 
ministically, the fastest convergence rate is obtained for 
t = 2/(En-i — Ei); semi- or fully stochastic multiplica- 
tions require smaller values of r to reduce the statistical 
noise @. 

Let x be an arbitrary initial vector satisfying 
(X {0) \i' {0) ) ^ 0. Then, repeated application of P to x (0) 
yields 



Px 



(0 



(2) 



and according to the power method, oc ip( ' for 

sufficiently large M. 
the expansion 



If coefficients wf^ are defined by 



N 



(3) 



the semistochastic representation of the weights wf* and 
the multiplication by P in Eqn. ||2J) are defined as follows. 

Let T> be the set of deterministic indices and S be 
the set of stochastic indices, where DU<S = {1, . . . , N}, 
V n S = 0, and \D\ <C N. Accordingly, P is the sum of 
a deterministic block P° , and a stochastic complement 
P 5 , 



P = P L 



P i 



where 



PP. 



V 



Pij, hj&'D, 
0, otherwise, 



(4) 



(5) 



and P s = P - P 

If the deterministic space is the entire space, then there 
is no sign problem or statistical noise. Consequently, it 
is expected that using a deterministic subspace, which 
is not the entire space, reduces the sign problem and 
statistical noise. 

The coefficients of the basis functions are represented 
as a population of walkers, where the number of walkers 
on an occupied \4>i) is 



rii = max(l, [KID, 



(6) 



where denotes the nearest integer and each walker has 
signed weight w,/n,-. 

Next, we proceed to the multiplication by P which 
evolves the coefficients from time t to time 4+1. 

• To account for the off-diagonal elements in P s . for 
each walker on \4>i), a move to \<pj) ^ \<f>i) is made 



with probability T^. 
tributes 



A single walker on \<j>i) con- 



0. 



to 



(t) 



otherwise 



(7) 



to the signed walker weight on \4>j). The choice 
of T determines the probability that particular off- 
diagonal moves are made. In this work, the near- 
uniform choice of Booth, Thorn and Alavi is used 
■ To control sign problems present in our exam- 
ples, we use the initiator idea [H, which we gener- 
alized in that the initiator threshold increases with 
the number of steps taken since the last visit to the 
deterministic space fioj . 

To account for the diagonal elements in P s , the 
contribution to the total signed walker weight on 
\4>j), with j G S, is 



P. 



j j 



w 



(8) 



• Deterministic evolution is performed with P^. The 
contribution to the signed weight on \<j>j), with j G 
V, is 



(9) 



P° is stored and applied as a sparse matrix. 



Finally, for each \4>j), all signed walker weight gen- 
erated on is summed, taking into account the 
sign of the contribution. To avoid the large compu- 
tational and memory cost of having small weights 
on a large number of basis states, basis states with 
weight less than some minimum cutoff . w m i n , are 



combined via an unbiased prescription [11 



After sufficiently many multiplications by P, contribu- 
tions from subdominant eigenvectors die out on average. 
At this point, the collection of averages begins. The most 
commonly employed estimator for the dominant eigen- 
value is the mixed estimator 



(V> (0) |#|V<T> 



(10) 



where the trial state \4>t) satisfies (i/A * 1 \iPt) / 0. 

The trial state \iPt) is a linear combination of basis 



states lli 



(ii) 



where T is the set of indices of those basis functions that 
contribute to the trial state. We require that |T| N, 
but not necessarily that TcP. 
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At any particular time t, the stochastic representation 
of the dominant eigenvector is 

i^°>>« J2 w fi&>> ( 12 ) 

where is the set of indices of basis functions oc- 

cupied by walkers at time t. The full representation of 
the dominant eigenvector is obtained by summing over 
Monte Carlo generations 

i^ (0) )-^E E «^i«fc>> as) 

gcn t=i ieww 

where N gen is the number of times P is applied after 
equilibration. 

For the trial state in Eqn. (fTTj) , 

p _ J2t=T Dig ww w z Sjgr Hjjdj , . 

Since -E m i X is a zero-variance-zero-bias estimator when 
I^t) is equal to the dominant eigenvector, improving the 
quality of \iPt) reduces fluctuations and bias in the mixed 
estimate of the dominant eigenvalue. This reduction can 
be achieved with almost no additional computational cost 
by storing nonzero X)jgT Hijdj terms. 

Applications. — The semistochastic method is now ap- 
plied to compute the ground state energy of both the 
simple-square 8x8 fermionic Hubbard model and the 
carbon dimcr. In both cases, we represent H in the basis 
of determinants formed from the restricted Hartree-Fock 
orbitals. For the former, these orbitals are the momen- 
tum eigenstates. For the latter, these orbitals are ob- 
tained by solving the restricted Hartree-Fock equations 
in cc-pVTZ basis set [l3|. The majority of the Hubbard 
calculations are performed for U/t = 4 where U is the 
on site Coulomb repulsion and t is the nearest neigh- 
bor hopping. This parameterization is considered to be 
in the intermediate coupling regime (the nonintcracting 
bandwidth being 8t), and has been used widely in the 
literature (l4j |. 

The trial wavefunction space and the deterministic 
space are generated with identical iterative schemes, but 
possibly different parameters. At each iteration, first de- 
fine a reference space as all states obtained in the pre- 
vious iteration. Second, generate a space which includes 
all determinants connected to the reference space by a 
single application of the Hamiltonian. Third, find the 
dominant eigenvector in this space. Fourth, truncate the 
space using a criterion based on the magnitude of the co- 
efficient of each state in the eigenvector. This truncated 
space becomes the reference for the next iteration. The 
reference for the first iteration is the Hartree-Fock state. 

For various sizes of the deterministic space, we demon- 
strate the improvements of SQMC over the purely 



stochastic method defined by a deterministic space which 
includes only the Hartree-Fock determinant. The purely 
stochastic method is almost the same as i-FCIQMC [7|, 
[H, aside from some details such as the use of real 
walker weights versus the integer walker weights used in 
FCIQMC and the use of a graduated initiator in SQMC. 
The most dramatic benefit of SQMC is in the efficiency, 
which is defined to be proportional to the inverse of the 
time required to obtain the ground state energy to a spec- 
ified level of uncertainty. 

Fig. [T] shows the efficiency gain of SQMC vs. the size 
of the deterministic space for the simple-square 8x8 
Hubbard model with U/t = 4 and 10 electrons. The 
orders of magnitude increases in efficiency demonstrate 
not only the benefits of SQMC, but also the benefits that 
can be achieved by improving the trial wavefunction. The 
efficiency benefit of just using the largest deterministic 
space is a factor of 22, while the benefit of just using the 
largest trial wavefunction is a factor of 28. Together the 
benefit is over 600 as seen in the plot, but the benefits 
need not be multiplicative for all systems. 
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FIG. 1. Efficiency of SQMC vs. dimension \T>\ of the de- 
terministic space for the simple-square 8x8 Hubbard model 
with U/t = 4 and 10 electrons. Results are shown for trial 
wavefunctions of increasing size. All values are normalized by 
the efficiency of the stochastic method (\T>\ — 1) with |T| = 1. 
The inset shows the |T| = 1 curve on an expanded scale. For 
this system, N w 10 12 . 

Fig. [2] shows the efficiency gain of SQMC vs. filling 
fraction for the simple-square 8x8 Hubbard model with 
U/t = 4. The efficiency gains increase with increasing 
filling fraction, demonstrating the potential for SQMC 
to study strongly correlated systems. 

SQMC produces large efficiency gains for chemical sys- 
tems as well. Fig.[3]shows the efficiency gain of SQMC vs. 
the size of the deterministic space for the carbon dimcr 
with a cc-pVTZ basis set [HI • The bottom two curves are 
for V and T generated with one applications of our itera- 
tive scheme which generate single and double excitations 
only. The largest efficiency gain for these is about 40. 
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FIG. 2. Ratio of the efficiency of SQMC to the stochastic 
method at the same filling vs. filling fraction for the simple- 
square 8x8 Hubbard model with U/t = 4. The trial wave- 
function for each of these calculations is the Hartree-Fock 
determinant. The deterministic space consists of states con- 
nected to the Hartree-Fock determinant by one application of 
the Hamiltonian, resulting in spaces of sizes 1412, 4088, 7424, 
14160, 16540, respectively. N ranges from roug hly 10 12 to 
10 35 . 



The top two curves are for V and T generated with two 
applications of our iterative scheme and, hence, include 
several chemically relevant quadruple excitations which 
are important for correctly describing the ground state 
wavefunction. The largest efficiency gain now jumps to 
over 1000! 
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FIG. 3. Efficiency of SQMC vs. dimension \T>\ of the deter- 
ministic space for the carbon dimer with a cc-pVTZ basis. 
Results are shown for trial wavefunctions of increasing size. 
The top two curves are for T> and T generated with two appli- 
cations of our iterative scheme and, hence, include chemically 
relevant quadruple excitations. All values are normalized by 
the efficiency of the stochastic method (\T>\ = 1) with \T\ = 1. 
For this system, N » 10 9 . 



Not only is SQMC much more efficient than the 
stochastic method, but in some cases, the initiator bias is 



significantly reduced. Fig.[4]shows the biased estimates of 
the energy as obtained by both the SQMC and stochastic 
method vs. the average number of occupied determinants 
for the 8x8 Hubbard model with U/t — 1 and 50 elec- 
trons. SQMC has essentially no bias. A larger average 
number of occupied determinants corresponds to using 
a larger walker population in the calculation. The time 
required for a calculation is proportional to the walker 
population. 
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FIG. 4. Energy of SQMC and the stochastic method vs. 
the average number of occupied determinants for the simple- 
square 8x8 Hubbard model with U/t = 1 and 50 electrons. 
The trial wavefunction for each of these calculations is the 
Hartree-Fock determinant. The deterministic space consists 
of states connected to the Hartree-Fock determinant by one 
application of the Hamiltonian, yielding a deterministic space 
of 16540 determinants. For this system, N « 10 35 . 

The reduction in initiator bias is not always large. 
Fig. [5] shows both the SQMC and stochastic method en- 
ergy vs. the average number of occupied determinants for 
the 8x8 Hubbard model with U/t = 4 and 10 electrons. 
SQMC has a reduced initiator bias for a small, but not 
for a large number of occupied determinants. However, 
for this system and all other systems studied SQMC has 
a smoother bias than the stochastic method. 

Conclusion. — The semistochastic power method, a 
hybrid with deterministic and stochastic components, 
was introduced for finding the dominant eigenvalue and 
sampling the corresponding eigenvector of a matrix. We 
showed that introducing a deterministic component into 
the stochastic method significantly reduces the statisti- 
cal noise while it does not compromise the ability of the 
method to deal with matrices well beyond the size that 
can be handled by deterministic methods. In particular, 
matrices ranging in order from 10 9 to 10 35 were tackled. 
Besides being more efficient than a purely stochastic ap- 
proach, the semistochastic method has in some cases the 
additional benefit of a much reduced initiator bias. The 
applications presented here are to systems with a sign 
problem, but the efficiency benefits of a semistochastic 
implementation of the power method extend to systems 
without a sign problem. 
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FIG. 5. Energy of SQMC and the stochastic method vs. 
the average number of occupied determinants for the simple- 
square 8x8 Hubbard model with U/t = 4 and 10 electrons. 
The trial wavefunction for each of these calculations is the 
Hartree-Fock determinant. The deterministic space reference 
state for each SQMC calculation is the Hartree-Fock deter- 
minant, yielding a deterministic space of 1412 determinants. 
For this system, N ^ 10 12 . 
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